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Abstract 



A class of second order approximations, called the weighted and shifted Griinwald difference op- 
erators, are proposed for Riemann-Liouville fractional derivatives, with their effective applications to 
numerically solving space fractional diffusion equations in one and two dimensions. The stability and 
' convergence of our difference schemes for space fractional diffusion equations with constant coefficients 

^ | . in one and two dimensions are theoretically established. Several numerical examples are implemented 

to testify the efficiency of the numerical schemes and confirm the convergence order, and the numerical 
results for variable coefficients problem are also presented. 

Keywords: Riemann-Liouville fractional derivative, Fractional diffusion equation, Weighted and 
shifted Griinwald difference operator. 

AMS subject classifications: 26A33, 65L12, 65L20 

>: 

^ ; 1 Introduction 

^ . Fractional calculus is a fundamentally mathematical tool for describing some special phenomenons arising 

from engineering and science [15, 18, 22]. One of its most important applications is to describe the subdif- 
fusion and superdiffusion process [5, 10, 16]. The suitable mathematical models are the diffusion equations 
with time and/or space fractional derivatives, where the classical first order derivative in time is replaced by 
the Caputo fractional derivative of order a € (0, 1), and the second order derivative in space is essentially 
replaced by the Riemann-Liouville fractional derivative of order a € (1,2]. The physical interpretation and 
practical applications of fractional diffusion equations have been discussed a lot with some common ideas 
[1, 9, 14]. Based on these, our main purpose of this paper is to study the higher accurate numerical solution 
of the space fractional diffusion equation by a novel finite difference approximation. 

From the perspective of the numerical analysis, there are some fundamental difficulties in numerically 
approximating the fractional derivatives, because some good properties of classical approximating opera- 
tors are lost. Over the last decades, the finite difference method has some developments in solving the 
fractional partial differential equations, e.g., [2, 12, 13, 27]. The Riemann-Liouville fractional derivative 
can be discretized by the standard Griinwald-Letnikov formula [18] with only the first order accuracy, but 
the difference scheme based on the Griinwald-Letnikov formula for time dependent problems is unstable 
[12]. To overcome this problem, Meerschaert and Tadjeran in [12] firstly proposed the shifted Griinwald- 
Letnikov formula to approximate fractional advection-dispersion flow equations. Recently, second order 
approximations to fractional derivatives are studied, Sousaa and Li presented a second order discretiza- 
tion for Riemann-Liouville fractional derivative and established an unconditionally stable weighted average 
finite difference method for one-dimensional fractional diffusion equation in [23], and the results in two- 
dimensional two-sided space fractional convection diffusion equation in finite domain can be seen in [6]; 
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Ortigueira [17] gave the "fractional centred derivative" to approximate the Riesz fractional derivative with 
second order accuracy, and this method was used by Celik and Duman in [2] to approximate fractional dif- 
fusion equation with the Riesz fractional derivative in a finite domain. In this paper, we propose a more 
general and flexible approach to approximate the Riemann-Liouville fractional derivative via combining the 
distinct shifted Grunwald-Letnikov formulae with their corresponding weights, and the weighted and shifted 
Griinwald-Letnikov formulae achieve second and higher order accuracy. A detailed algorithm shows that 
the weights are related to not only the shifted numbers but also the order of the fractional derivative, which 
implies the numerical algorithm is more related to the equation itself. 

The paper is briefly summarized as follows. In Sec. 2, we propose a class of discrete operators to 
approximate the Riemann-Liouville fractional derivatives with high order truncating errors. In Sec. 3 and 
4, one dimensional and two dimensional fractional diffusion equations are numerically solved by using the 
finite difference method based on the weighted and shifted Grunwald-Letnikov formulae, and the stability 
analysis of each case is presented. We prove that the finite difference solutions approximate the exact ones 
with 0(r 2 + h 2 ) in the discrete L 2 norm. Some numerical experiments are performed in Sec. 5 to verify the 
efficiency and accuracy of the methods. And the concluding remarks are given in the last Section. 

2 High Order Approximations for Riemann-Liouville Fractional Derivatives 

We begin with the definitions of the Riemann-Liouville fractional derivatives and the properties of their 
Fourier transform. 

Definition 1 ([18]). The a(n — 1 < a < n) order left and right Riemann-Liouville fractional derivatives of 
the function u{x) on [a, b] are defined as 

(I) left Riemann-Liouville fractional derivative: 



If a = n, then a D%u(x) = £^u(x) and x D^u(x) = (-l) n ^u(x). 

Property 1 ([8]). Let a > 0, u € C^°(0), Slci The Fourier transforms of the left and right Riemann- 
Liouville fractional derivatives satisfy 




(2) right Riemann-Liouville fractional derivative: 




where u(ca) denotes the Fourier transform ofu, 




In [12], the shifted Grunwald difference operator 




(2.1) 



fc=0 
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approximates the Riemann-Liouville fractional derivative uniformly with first order accuracy, i.e., 

A%pu{x) = -ooD^u(x) + 0(h), (2.2) 

where p is an integer and = (— l) fc (^). In fact, the coefficients in (2.1) are the coefficients of the 
power series of the function (1 — z) a , 



(i-*r=£(-i)*ry=£*£ 

fc=0 V 7 fc=0 

for all |z| < 1, and they can be evaluated recursively 



(a) 



(a) (-. a+l\ ( a ) 

56 5fc - (1 — )9k-v k = i' 2 --- 

Lemma 1 ([12, 13, 18]). 77ie coefficients in (2.1) satisfy the following properties for 1 < a <2, 
( = 1, = -a < 0, 

oo m 

E^ a) = ' £<4 a) <o,m>i. 

k fe=0 k=0 



(2.3) 



(2.4) 



(2.5) 



2.1 Second Order Approximations 

Inspired by the shifted Griinwald difference operator (2.1) and multi-step method, we derive the following 
second order approximation for the Riemann-Liouville fractional derivatives. 

Theorem 1. Let u £ L 1 (R), _ooD° +2 ti and its Fourier transform belong to L 1 (R), and define the weighted 
and shifted Griinwald difference (WSGD) operator by 



then we have 

L V^ q u(x) = _ QO D«u(x) + 0(h 2 ) 
uniformly for x € R, where p, q are integers and p ^ q. 

Note. The role of p and q is symmetric, i.e., L^hpq u ( x ) = q P u ( x )- 
Proof of Theorem 1. By the definition of A^ p in (2.1), we can rewrite the WSGD operator as 

a -2q 1 ^ {a) 2p-al^ {a) 



(2.6) 



(2.7) 



L^h,p,qU(x) 



2(p - q) h< 



Y / g^u(x-(k-p)h) + 



k=0 



2{p - q) h° 



Y^gr^x-ik-q)^. (2.8) 



k=0 



Taking Fourier transform on (2.8), we obtain 

-. oo „ 
& [ L V h ^ q u\ {^-—^ g k [ 2{p _ fl ^ 



h c 

k=0 

1 / a - 2q 



2p-a Mk _ q)h 
2(p-q) 



u(uj) 



h«\2(p-q) 



(2.9) 



'{ ^ — \w p {iojh) + ^ ° W q (iujh) )u(cj), 



2(p - q) 



2(p - q) 



where 

W r (z) = (^—y e rz = l + (r-^)z + 0(z 2 ), r = p,q. (2.10) 
Denoting $(u>, h) = &[L^ pq u\{u) - ^[-ooD^u]^), then from (2.9) and (2.10) there exists 

\4>(u,h)\ < Ch 2 \iu\ a+2 \u{uj)\. (2.11) 
With the condition ^[_ OOJ D^ +2 u](a;) G it yields 

IlZ^u - -oc^>| = |0| < ^ | < C||^[_oo££ +2 u]MI|li& 2 = 0(/i 2 )- (2-12) 

□ 

Remark 1. For the right Riemann-Liouville fractional derivative, similar to Theorem 1, we can check that 
R Vl M u(x) = ^-^B^ p u(x) + ^^BlJ{x) = x D^u(x) + 0(h\ (2.13) 

uniformly for x € M under the conditions that u € L 1 (IR), x D^ 2 u and its Fourier transform belong to 
L 1 (M), where p, q are integers and 

1 oo 

B hA x ) = ^ + r)h). (2.14) 

k=0 

Remark 2. Considering a well defined function u{x) on the bounded interval [a, b], ifu(a) = or u(b) = 0, 
the function u(x) can be zero extended for x < a or x > b. And then the a order left and right Riemann- 
Liouville fractional derivatives of u(x) at each point x can be approximated by the WSGD operators with 
second order accuracy 

A m+p A [^i+a 

a D«u{x) = ^ 9 { k a) n(x-(k-p)h) + 1 ^ Yl 9 { k a) n(x-(k-q)h)+0(h 2 ), 

A i b -irl+P A lnr\+i 

x DZu(x) = 1 ± Yl 9 k a) u(x+{k-p)h) + ^ Y, 9 k a) n(x+(k-q)h)+0(h 2 ), 

k=0 k=0 

where Ai = g^fy, A 2 = |^y. 

Remark 3. The integers p, q are the numbers of the points located on the right/left hand of the point x used 
for evaluating the a order left/right Riemann-Liouville fractional derivatives at x, thus, when employing the 
difference method with (2.15) for approximating non-periodic fractional differential equations on bounded 
interval, p, q should be chosen satisfying \p\ < l,\q\ < 1 to ensure that the nodes at which the values of u 
needed in (2.15) are within the bounded interval; otherwise, we need to use another way to discretize the 
fractional derivative when x is close to the right/left boundary. When (p, q) = (0, —1), the approximation 
method turns out to be unstable for time dependent problems. So two sets of (p, q) can be selected to estab- 
lish the difference scheme for fractional diffusion equations, that is (1,0), (1, —1), and the corresponding 
weights in (2.6) and (2.13) are (§, 2 =f) and (^±2, 2=2). For a = 2, the WSGD operator (2.6) is the 
centered difference approximation of second order derivative when (p,q) equals to (1,0) or (1,-1); for 
a = 1, (p,q) = (1, 0), the centered difference scheme for first order derivative is recovered. 
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The simplified forms of the discreted approximations (2. 15) for Riemann-Liouville fractional derivatives 
with (p,q) = (1,0), (1,-1) are 

j i+l 

a D*u(xi) = j^Yl w k ] ui x i-k+i) + 0{h 2 ), 

k=0 

x D%u(xi) = — ^2 w l u ( x i+k-i) + 0(h 2 ), 



(2.16) 



fc=0 



where 



fog) = (1,0), < 



(n) 



O: 



(«) ...(«) 



n 



lo^ -I- - — n w A- > 1 • 



2 - a 



,(«) 



<M)-(1,-1). 4°' = ^""US"' = 



(2.17) 



With Lemma 1 and some calculations, we obtain the properties of the coefficients w^* in (2.16) correspond- 
ing to (p, q) = (1, 0), (1, —1) as follows. 

Lemma 2. The coefficients in (2.16) satisfy the following properties for 1 < a < 2, 
(1) if(p,q) = (1,0), 

i a ) a f a ) 2-a-a 2 r a ) a(a 2 + a - 4) 

= o ' ™1 = o < ' ^2 = A ' 



J ° " 2' 1 " 2 
1 > > wi a) >w^>...>0, 

oo m 

J2w^=0, 2>«<0,m>2; 



(2.18) 



fc=0 fc=0 

(2) (f(p, g ) = (1,-1), 



(a) 2 + a ( a ) 
= ; , W) 



2a + a 2 



<0, 



( a ) a 3 + a 2 — 4a + 4 ( a ) a(2 — a)(a 2 + a 

W \ = >o, w y = 

1 > > w { 2 a) > w'f > wi? ] > . . . > 0, 

oo m 

= 0, < 0, m = 1 or m > 3. 



< 0, 



(2.19) 



^ fc=0 



k=0 



Next, we will explore the properties of the eigenvalues of the difference matrix of (2.16) on grid points 
{x^ = a + kh, h = (b — a)/n,k = 1, 2, . . . ,n — 1}. In the following, we denote by i7 the symmetric 
(respectively, hermitian) part of A if A is real (respectively, complex) matrix. 

Lemma 3 ([20]). A real matrix A of order n is positive definite if and only if its symmetric part H = A+ 2 A 
is positive definite; H is positive definite if and only if the eigenvalues of H are positive. 

Lemma 4 ([20]). If A e C nxn , let H = A+ 2 A * be the hermitian part of A, A* the conjugate transpose of A, 
then for any eigenvalue X of A, there exists 

Knn(H) < Re (A) < A max (#), 

where Re(A) represents the real part of X, and X m i Q (H), X max (H) are the minimum and maximum of the 
eigenvalues of H. 
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Definition 2 ([4]). Let Toeplitz matrix T n be of the following form, 



( h t_i • • • £2-? 

tl to t-i 

: ti to 

tn-2 

\t n -l tn-2 ■ ■ ■ h 
n-1 



h-n\ 
t-1 

to ) 



if the diagonals {tfc}^ = „ n+1 are the Fourier coefficients of a function f, i.e., 
then the function f is called the generating function ofT n . 

Lemma 5 (Grenander-Szego theorem [4, 3]). For the above Toeplitz matrix T n , if f is a 2ir-periodic con- 
tinuous real-valued function defined on [—it, it], denote A mm (T n ) and A max (T n ) as the smallest and largest 
eigenvalues ofT n , respectively. Then we have 

fmm — A m ; n (T' rl ) < A max (T' rl ) < /maxi 

where fmm, /max denote the minimum and maximum values of f(x). Moreover, if / m ; n < / max » then all 
eigenvalues ofT n satisfy 

fmm ^ A(X' n ) < /max; 

for all n > 0; and furthermore if fmm > 0, then T n is positive definite. 
Theorem 2. Let matrix A be of the following form, 



A 



?<>, 



(a) 



W 



(a) 


(«) ..,(«) 



w 







w. 



(«) ..,(«) 



(a) 



n-2 



(a) 



(a) 
n-2 



W 



\ 



(a) 

(a) (a) . 
W l / 



(2.20) 



where the diagonals {w^}^_q are the coefficients given in (2.16) corresponding to (p,q) = (1,0) or 
(1, —1). Then we have that any eigenvalue A of A satisfies 

(1) Re(X) = 0, for (p,q) = (1,0), a = 1, 

(2) Re (A) < 0,for(p,q) = (1,0), 1 < a < 2, 

(3) Re (A) < 0, /or (p, 9 ) = (1, -1), 1 < a < 2. 

Moreover, when 1 < a < 2, matrix A is negative definite, and the real parts of the eigenvalues of matrix 
c\A + C2A 1 - are less than 0, where c±,C2 > 0, c\ + c\ ^ 0. 

Proof. We consider the symmetric part of matrix A, denoted as H = A+ 2 AT ■ The generating functions of A 
and A T are 



/a(x) = £ 



(«) i(fc-l)a 



(a) -i(k-l)x 



k=0 



k=0 
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respectively. Then f(a;x) = ^j_E_j^Al_g_ j s generating function of H, and f(a;x) is a periodic 
continuous real- valued function on [— tt, it] since /a(x) and /^t(x) are mutually conjugated. 
Case (p, q) = (1, 0): with the corresponding coefficients given by (2.17), then 



k=0 k=0 

. oo „ oo oo „ oo 



2 V 9 / J fe ~ ' 2 / ^ J fc ~ 2~ / ^ J fc ~ 1 2 ' ' 

fc=0 fc=0 fc=0 fc=0 

= I (V J:r (l - e iX ) a + e lx {l - e~ lx ) a ^j + ((1 - e ix ) Q + (1 - e" 

Next we check f(a;x) < for 1 < a < 2. Since f(a;x) is a real-valued and even function, we just 
consider its principal value on [0, tt]. By the formula 



e w _ e «P = 2% sin ( — 2 , 
we obtain 

/(a; x) = (2 sin(|)) a (| cos - tt) - s) + ^ cos (|(x - vr))) . (2.21) 

It is easy to prove that f(a; x) decreases with respect to a, then f(a; x) < /(l; x) = 0; by Lemma 4 and 5, 
Re(A) = for a = 1, and /(a; x) is not identically zero for 1 < a < 2, then we get Re(A) < 0. 

Case (p, q) = (1, —1): the corresponding generating function f(a; x) of A+A can be calculated in the 
following form with coefficients wf 1 given by (2.17), 



1 oo oo 



(a) -i(k-l)x 



w k 'e 

k=0 k=0 

oo oo _ oo oo 

L-** g { ^e ikx + e ix ^ gjpe**) + (f x E Sfc*^** + e ~" E ^e"** 

fc=0 k=0 k=0 k=0 

* r e~ ix (l - e ix ) a + e ix (l - e~ ix ) a ) + ^—^ (e ix (l - e ix ) a + e~ ix (l - e~ ix ) a 



Next we check f(a;x) < for 1 < a < 2. Since f(a;x) is a real-valued and even function, we just 
consider its principal value on [0, tt]. By simple calculation, we obtain 

f(a;x) = (2sin(-)) l^- sin (-(x - vr)) sin(x) + cos (-(x - tt)) cos(x)J . (2.22) 

We can also check that f(a; x) decreases with respect to a, then f(a; x) < /(l; x) = —2 sin 4 (|) < 0, then 
by Lemma 4 and 5, we get Re(A) < for 1 < a < 2. 

From the above discussions and Lemma 5, we know, for 1 < a < 2, the matrix \{A + A T ) is negative 
definite, which implies matrix A is negative definite by Lemma 3. And the symmetric part of matrix c\A + 
c 2 A T is ^4fi(A + A T ), thus we obtain Re(A(ciA + c 2 A T )) < for 1 < a < 2. □ 

Remark 4. For the case (p,q) = (1, 0) awe? 1 < a < 2, we can check that the symmetric part H of matrix 
A in (2.20) is strictly diagonally dominant by using Lemma 2, and the elements of the main diagonal of H 
are negative, then the eigenvalues of H are less than zero by the Gershgorin circle theorem ([20],P188), 
therefore, with Lemma 3 and 4, we can also get Re(A(A)) < 0, and A is negative definite. 

Remark 5. By the same approach described in Theorem 2, we can verify that the generating function of the 
symmetric part of difference matrix for (p, q) = (0, —1) is not identically negative when 1 < a < 2, which 
leads to the instability of the difference method to fractional diffusion equations for the same reason in the 
stability analysis in the sequel. 
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2.2 Third Order Approximations 

Similar to the second order approximations for Riemann-Liouville fractional derivatives, we give a combi- 
nation of three shifted Griinwald difference operators 

LQt PA , r u{x) = XiAl p u{x) + X 2 Al q u{x) + X 3 Al r u(x), (2.23) 

where p, q, r are integers and mutually non-equal, and 

Ylqr - (6q + 6r + l)a + 3a 2 



A 



12(qr — pq — pr + p 2 ) 

3o. 

I2pq - (6p + 6q + l)a + 3a 2 



\2pr- (6p + 6r + l)a + 3a 2 

A2 = To? ; 2\ ' (2.24) 

\l\pr — pq — qr + g z j 



3 12(pg — pr — qr + r 2 ) 

Assuming u £ and taking Fourier transform on (2.23), we get 

■^b^,,,r"](w) = (iw) a (A 1 VF p (^/ l ) + A 2 W 3 (ia;/t) + X 3 W r {iojh)ju(co) 
= {iuj) a (l + C(iujh) 3 ^u(uj), 

where W s (z) is defined in (2.10). If _ooL)" +3 w and its Fourier transform belong to L 1 (R), then we have 



(2.25) 



\Lyh,p,q,r u 



,D°u\ < ±- [ \&[ L Gt m , r u - -ooD«u}\ 

Z7r Jr {I. lb) 

< Cll^^oo^nlHII^/i 3 = 0{h 3 ). 



The above results can be stated in the following theorem. 

Theorem 3. Let u € L 1 (M), -^D^^u and its Fourier transform belong to and the following 

3-WSGD operator (2.23) satisfies 

LQh, P , q ,M x ) = -ocD%u(x) + 0(h 3 ), (2.27) 
uniformly for x € R. 

If u £ L X (R), x D^~ 3 u and its Fourier transform belong to L 1 (R), we also have 

uniformly for x G R, where the operator Bff s is given by (2.14), and Aj, i = 1, 2, 3 are the same as (2.24). 

As stated in Remark 3, the 3-WSGD operator can be utilized for approximating Riemann-Liouville 
fractional differential equations on bounded domain by finite difference method when choosing (p, q, r) = 
(1,0, —1), then the corresponding weight coefficients in (2.24) are Ai = J^a + |a 2 , A2 = 1 + j^a — 
\o? , A3 = — + \ot 2 . For function u(x) satisfying u(a) = u(b) = on grid points {x^ = a + kh, h = 
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(b — a)/n, k = 1, 



1}, the approximation matrix of (2.23) with (p, q, r) = (1,0, —1) is 



G=Ai 



/ (") 
' 9i 

(a) 
52 



(a) 

\a {a) 

\9 n -l 



(a) 
50 
(a) 

9i 

(a) 
92 



(a) 
(a) 

9i 



(a) 
9n-2 



(a) 
92 



\ 



(a) 

So 

(a) , 

5i 7 



+ A 2 



\ 



(«) 



(a) 
\9n-2 



(a) 
50 



(a) 



.9i 



(a) 



(a) 
#0 



(«) 



(a) , 

ffo 7 



+ A 3 



/ 

(a) 

9o 



(a) 
Sn-4 

\fln-3 



\ 



(2.29) 



(a) 

So 



(a) 
9n-A 



(a) 

So 



0/ 



Example 1. We utilize the approximation (2.23) for simulating the steady state fractional diffusion problem 



- D%u(x) 



T(3 + a) 



x e (0,1) 



(2.30) 



r 2+a 



with u(0) = 0, u(l) = 1, and 1 < a < 2. 77ie exac? solution is u(x) 

The 3-WSGD operator with (p, g, r) = (1, 0, —1) is utilized for computing the solution of Example 1, 
the numerical results are given in Table 1, from which the order and accuracy of the 3-WSGD operator is 



verified. As we do in the above, the generating function of the symmetric part 



G+G 1 



of the Toeplitz matrix 



Table 1: The maximum and L 2 errors and their convergence rates to Example 1 approximated by the 3- 
WSGD operator for a = 1.1, 1.9. 



iV 




a 


1.1 






a 


1.9 




||u n -*7"||oo 


rate 


\\u n - U n \\ 


rate 




rate 


\\u n - U n \\ 


rate 


8 


9.48629E-04 




5.92003E-04 




3.20333E-04 




1.59788E-04 




16 


1.19530E-04 


2.99 


7.51799E-05 


2.98 


2.29262E-05 


3.80 


1.04858E-05 


3.93 


32 


1.50130E-05 


2.99 


9.47995E-06 


2.99 


1.58500E-06 


3.85 


6.71546E-07 


3.96 


64 


1.88094E-06 


3.00 


1.18999E-06 


2.99 


1.07818E-07 


3.88 


4.24776E-08 


3.98 


128 


2.35382E-07 


3.00 


1.49052E-07 


3.00 


7.27733E-09 


3.89 


2.67067E-09 


3.99 


256 


2.94392E-08 


3.00 


1.86501E-08 


3.00 


4.89318E-10 


3.89 


1.67325E-10 


4.00 



Gis 

f(a; x) =(A a + l a 2) (e"«(l - e ix ) a + e fa (l - e 

+ ^ + ^ a_ r 2 K (1 - eiT+(1 - e_ 

+ (-^a + la 2 )(e-(l-e i T + e- 

x G [— 7T, it]. As matrix G+ ^ T is symmetric, thus f(a; x) is a real-valued and even function, so we consider 
it on [0, 7r] and get 

f(a;x) =(2sin(|)) Q ((^« + T^« 2 ) cos (f ( x ~ 7r ) ~ x ) + [\ + ^« ~ £« 2 ) cos (^(x - vr)) 
+ ( " ^« + Yq® 2 ) cos (f ( x " ^ + x )) • 



9 



We can check that f(a; x) is not identically positive or negative for 1 < a < 2, and the real parts of the 
eigenvalues of matrix G are not always negative, so the finite difference scheme using (2.23) or (2.28) for 
time dependent fractional problems will not be unconditionally stable. 

3 One Dimensional Space Fractional Diffusion Equation 

In this section, we consider the following two-sided one dimensional space fractional diffusion equation 



where both n D" and X D^ are Riemann-Liouville fractional operators with 1 < a < 2. The diffusion 
coefficients K\ and K 2 are nonnegative constants with K\ + K 2 7^ 0. And if K\ ^ 0, then 4> a {t) = 0; if 
K 2 7^ 0, then <pb(t) = 0. Next we will discretize the problem (3.1) by the second order accurate WSGD 
formula (2.16). In the analysis of the numerical method that follows, we assume that (3.1) has a unique and 
sufficiently smooth solution. 

3.1 CN-WSGD scheme 

We partition the interval [a, b] into a uniform mesh with the space step h = (b — a)/N and the time step 
r = T/M, where N, M being two positive integers. And the set of grid points are denoted by xi = ih and 
t n = nr for 1 < i < N and < n < M. Let t n+1 / 2 = {t n + i n+ i)/2 for < n < M — 1, and we use the 
following notations 




(3.1) 



< = u(xi,t n ), /• 



n+1/2 



n+l 



O/T. 



Using the Crank-Nicolson technique for the time discretization of (3.1) leads to 



S t uf - - (K x { a D«u)^ + ^( aJ D»™ +1 + K 2 ( x B?u)? + K 2 ( X D^ +1 ) = /f +1/2 + 0(r 2 ). 



In space discretization, we choose the WSGD operators lV^ p q u(x,t) and rD^ q u(x , t) to approximate 
the Riemann-Liouville fractional derivatives a D^u(x, t) and x D£u(x, t) with second order accuracy, re- 
spectively, and (p, q) = (1, 0) or (1, —1). This implies that 



,n+l/2 n 



where 



sf\<c(r 2 + h 2 ). 



(3.3) 



Multiplying (3.2) by r and separating the time layers, we have 




n 1 "A' <7-,cr n , "z 1 -r-.a „,n , ^ f 7 



.n+1/2 



+ 0(T 3 + rh 2 ). 



(3.4) 
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Substituting L^p^u, rD^^u by (2.16), we obtain that 



Tj r i+1 N-i+1 

u n+i _ fvt£ »„(«) a ,w+i _ K 2T ^ 



2h c 



k=0 
i+1 



k=0 
N-i+1 



fc=0 k=0 

Denoting U™ as the numerical approximation of uf, we derive the CN-WSGD scheme for (3.1) 

V i+1 XT N-i+1 

n+ l_AlTv- ( Q ) n+ l _K2T w ^)jjn+l 
i Oh" k U i-k+l nua k U i+k-1 



2h a ^ k 2h c 

k=0 

rrn , KlT ST ln ( a hr + K2T ly \^\n {a) n n + T f n+1/2 



fc=0 
N-i+1 



k=0 k=0 

For the convenience of implementation, using the matrix form of the grid functions 



/ j.n+l/2 j.n+1/2 ,.n+l/2\ T 
I t>l > ^2 ) - - - > U N-1 I > * = [Jl ifo >'■■ ,JN-1 ) ' 



makes the finite difference scheme (3.6) be described as 

1 " ih^ {KlA + K2AT) ) un+1 = ( I+ i {KlA + K2AT) ) un + rpn + RT 



where A is given by (2.20) and 



2h a 



K lW { 2 a) + K 2 wi a) ' 



(uz + ut 1 ) + 



2h° 



(a) 



K 2 wl 
lK lW { a) + K 2 wi a \ 



3.2 Stability and Convergence 

Now we consider the stability and convergence analysis for the CN-WSGD scheme (3.7). Define 

Vh = {v : v = {vi} is a grid function in {xi = ih} 1 ^ 1 and vq = vn = 0}. 

For any v = {v^} € Vh, we define its pointwise maximum norm 

Halloo = max \vA 
l<i<N-\ 

and the following discrete norm 

IM 



N-l 



(3.5) 



(3.6) 



(3.7) 



{U n N + U n N +1 ). (3.8) 



(3.9) 



Theorem 4. The finite difference scheme (3.6) is unconditionally stable. 
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Proof. Denoting B = -^{K\A + K2A T ). The matrix form of the difference approximation for problem 
(3.1) can be rewritten as 

(I - B)U n+1 = (I + B)U n + rF n + H n . (3.10) 

If denote A as an eigenvalue of matrix B, then is the eigenvalue of matrix (I — B)^ 1 ^ + B). The 
result of Theorem 2 shows that the eigenvalues of matrix B+ 2 gT = ^^^-(A + A T ) are negative, thus 
Re(A) < 0, which implies that |^r| < 1. Therefore, the spectral radius of matrix (I — + B) is less 

than one, then the discreted scheme (3.6) is unconditionally stable. □ 

Remark 6. Considering the 6 weighted scheme for the time discretization of (3.1), then the iterative matrix 
of the full discrete scheme is 

(l-9B)~\l + {l-6)B), (3.11) 

ifX is an eigenvalue of matrix B, then the eigenvalue of (3.11) is ~^fenr^ - As Re(X) < 0, it is easy to check 
that 



1 + (1-0)A 



< 1 (3.12) 



1-0A 

for g < 9 < 1. Then the 9 weighted WSGD scheme for (3.1) is unconditionally stable when ~ < 9 <1. 

Before verifying the unconditional convergence of the scheme (3.6), we first present the discrete Gron- 
wall's inequality. 

Lemma 6 ([19]). Assume that {k n } and {p n } are nonnegative sequences, and the sequence {4> n } satisfies 

n— 1 n— 1 

0o < 9o, <t>n < 9o + ^Pi + X! kl< ^ u n ~ 1 ' 

1=0 1=0 
where go > 0. Then the sequence {<fi n } satisfies 

n—l n—1 



(t>n<[go + ^2pi)exv{%2ki), n>l. (3.13) 

1=0 1=0 

Theorem 5. Let uf be the exact solution of problem (3.1), and U™ the solution of the finite difference scheme 
(3.6), then for all 1 < n < M, we have 

|K - U n \\ < c(t 2 + h 2 ), (3.14) 

where c denotes a positive constant and \\ ■ \\ stands for the discrete L 2 -norm. 

Proof Let e™ = — f7J\ and from (3.5) and (3.6) we have 



(e n+i - e n ) - -^-A(e n+L + e n ) - — —A' (e n+i + e n ) = re n , (3.15) 

where 

\, n — TT n n n — TT n ■ ■ ■ n n — TT n \ T c n — ( r n c n . . . r n \ 
"1 u l ) "2 u 2i > a N-l u N—l) ' e — \ 1 ' 2 ' ' N-l J 



T 
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Multiplying (3.15) by h, and acting (e n+1 + e n ) T on both sides, we obtain that 

h(e n+1 + e n ) T I(e n+1 - e n ) - J^(e n+1 + e n ) T A(e n+1 + e n ) 



2h a ~ 



2h a ~ l K 

By Theorem 2, A and its transpose A T both being the negative definite matrices, we get 

( e "+i + e ") T A(e n+1 + e «) < 0, (e n+1 + e n ) T ^ T (e n+1 + e n ) < 0. (3.16) 

And it yields that 

N-l N-l 

h E (( e " +1 ) 2 " ( e ") 2 ) < rh E + e?)e?. (3.17) 

i=l i=i 

Summing up for all < k < n — 1, we have 

N-l N-ln-1 JV-ln-1 iV-1 

^ E ( e *) 2 * Th E E ( e * +1 + e e * = ^ E E e * C^ 1 + ^ + ^ E ^r 1 

i=l i=l k=0 i=l k=l i=l .n\ 

, JV-ln-1 , N-l n-l „ 7 iV — 1 , N—l 

< y E E( e ') 2 + y E E + 4) + ^ E ( e ?) 2 + £ E (-r 1 ) 2 - 

i=l k=l i=l k=l i=l i=l 

By noting that ef < c(t 2 + h 2 ), and utilizing the discrete Gronwall's inequality, we obtain that 

n—l „ 2 

||e"|| 2 <rEl|e fc || 2 +(£(r 2 + / i 2 )) < exp(T) (c(r 2 + h 2 )) < c ((r 2 + h 2 ) 2 ) , (3.19) 
fc=i 

which is the result that we need. □ 

4 Two Dimensional Space Fractional Diffusion Equation 

We next consider the following two-sided space fractional diffusion equation in two dimensions 

+ (k^cD$u(x, y, t) + K^ y D p d u(x, y, t)) + f(x, y, t), (x, y, t) € Q x [0, T], 
u{x,y,0) =u {x,y), (x,y)eU, 
u(x,y,t) = ip(x,y,t), (x,y,t) Gffix [0,T], 

where = (a,b) x (c,d), a D^, x D^ and c Dy, y D^ are Riemann-Liouville fractional operators with 1 < 
a, (3 < 2. The diffusion coefficients satisfy K+, > 0, i = 1,2, (iq + ) 2 + (iv~ 2 + ) 2 ^ and (Iff) 2 + 
(iff) 2 7^ 0. And the boundary function tp satisfies, if iff / 0, then ip(a,y,t) = 0; if iff ^ 0, then 
ip(b, y, t) = 0; if iff 7^ 0, then tp(x, c, t) = 0; if iff ^ 0, then ip(x, d, t) = 0. We assume that (4.1) has a 
unique and sufficiently smooth solution. 
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4.1 CN-WSGD scheme 

Now we establish the Crank-Nicolson difference scheme by using WSGD formula (2.16) for problem (4.1). 
We partition the domain Q into a uniform mesh with the space steps h x = (b — a) /N x ,h y = (d — c) /N y and 
the time step r = T/M, where N x , N y ,M being positive integers. And the set of grid points are denoted by 
Xi = ih x , yj = jh y and t n = nr for 1 < i < N x , 1 < j < N y and < n < M. Let t n+1/2 = (t n +t n+1 )/2 
for < n < M — 1, and we use the following notations 

u i,j = u {xi,yj,t n ), f™j 1/2 = f{xi,yj,t n+l/2 ), Stu^j = (u^ 1 - u\j)jT. 
Discretizing (4.1) in time direction leads to 

5 t ul = I (Kt( a D a x u)^ + K+( x DZu)^ + K-UD^ 1 + K^^u)^ ^ 
+ K+iaD-u)^ + K+{ x Dtu)l 3 + KTUDjtOS,- + K^( y D^) + f^ 1 ' 2 + 0(r 2 ). 



In space discretization, we choose the WSGD operators LVh x p <? n ' R^h x v i u anc * L h y p q u ^ R h y p q u t0 
respectively approximate the fractional diffusion terms a D x u, x D^u and c D y u, y D^u. And multiplying 
(4.2) by r and separating the time layers, we have that 



K i T -ry* K^r <r><* K i T t>P K * T T>P 

2 LU h x ,p,q 2 R h ^P,1 2 h V'P'1 2 /"// '•' 

1 + ^-i P ^,P,g + -^-R V h x ,p,q + + —^-R V hy !P! q) U i,j + T 4j + h3 ' 

(4.3) 

where |e™ ■ | < c(r 2 + /i 2 ) denotes the truncation error. And we denote 

£ = K+ L Vl )M + K$ R Vt M § = K^ L vi ypq + K 2R vl ypq . 
For simplicity, the step sizes are chosen as the same, h = h x = h y . Using the Taylor expansion, we have 



t 2 „, _ r 3 



-T-^/K;' »;:,; = - I {(K+ a D^ + K+ x D^(K^ c D^ + K^ y D^ Ut ) ij + 0(r 5 + T 3 h 2 ). (4.4) 
Adding formula (4.4) to the right-hand side of (4.3) and making the factorization leads to 



(l " = (l + (l + T ^)< 3 + rf^ 1/2 + rel 3 + 0(r 3 + r 3 /, 2 ). (4.5) 

Denoting by U™j the numerical approximation to u™j, we obtain the finite difference approximation for 
problem (4.1) 

1 " \S£) (l - 1^)U^ = (l + I#) (l + 1^)^. + r^ 2 . (4.6) 
For efficiently solving (4.6), the following techniques can be used. Peaceman-Rachford ADI scheme [24]: 



1 " 2<)^ = (l + 2^)^ + 2^ ' (4 " 7a) 

> - i^)^ 1 = (i + + W 1 ' 2 - (4 - 7b) 
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Douglas ADI scheme [7]: 



1 " M V% = ( 1 + -5 a x + r5 J) + r/^ 1/2 , (4.8a) 



2 x ) *J V 2 

l-^)^^^-^^- (4.8b) 
D'Yakonov ADI scheme [24]: 

(l " ^) ^ = (l + \?>i) (l + ^) + r/^ 1/2 , (4.9a) 



!- ^J^/ (4 - 9b) 
A simple calculation shows that 

^SStfjft 1 ' 3 = ^{KtaD a x + Kt x DZ){K^ c Dl + K^D^f^ 2 + 0(r 3 fc 2 ). (4.10) 
Then from (4.5) and (4.10), it yields that 

(l " \%) (l " = (l + (l + ^)<, + r/^ 1/2 + + re?,. (4-11) 
where 

^ = e?j " ^(K+ a D% + K+ x D?)(Kt c D$ + K~ y D P d )fZ+ 1/2 + 0(r 2 + r 2 ^ 2 ). (4.12) 
Eliminating the truncating error and denoting [/"• as the numerical approximation of u" -, we have 

( x - iO ( x - i^)^ 1 = ( x + ( x + + <? 1/2 + T^tt 1 ' 2 - (4 - 13) 

Introducing the intermediate variable V^™ , we obtain the locally one-dimensional (LOD) scheme mentioned 
in [21, 28], 

(l - T -5-)v- = (l + + I(l + jjf' 2 , (4.14a) 

(i - = (i + \4)v% + K 1 - ^W 1/2 - < 4 - i4b) 

4.2 Stability and Convergence 

Now we consider the stability and convergence analysis for the CN-WSGD scheme (4.6). Define the sets of 
the index of the interior and boundary mesh grid points in domain [a, b] x [c, d], respectively, as 

Ah = : 1 < i < N x - 1, 1 < j < N y - 1}, 

8A h = : i = 0,N X ;0< j < N y } U :0<i<N x ;j = 0,N y }. 

For any v = {v{\ € Vh, we define its pointwise maximum norm and discrete L 2 norm, respectively, as 
follows 



N X ~lNy-l 

v E E v h ( 4 - 15 ) 

i=l j=l 



^ \\oo liAU, - rt - r-^,?? \ 
(iJ)€A h ^ 

where 

Vh = {v : v = {vij} is a grid function in and Vij = on dAh}- 
In the following, we list some properties of Kronecker products of matrices. 
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Lemma 7 ([11]). Let A G R nxn have eigenvalues {Xi}f =1 , and B G R mxm have eigenvalues {fj,j}^L v 
Then the ran eigenvalues of A® B, which represents the kronecker product of matrix A and B, are 

Lemma 8 ([11]). Let A G R mxn ,B G R rxs ,C G R nxp ,D G R sxt . Then 

(A ®B)(C®D) = AC® BD (g R mrxpt ). (4.16) 
Moreover, if A, B G R nxn , I is a unit matrix of order n, then matrices I <g) A and B <g> I commute. 
Lemma 9 ([1 1]). For all A and B, (A ® B) T = A T ® 5 T . 

In the theoretical analysis of the numerical method, we choose N = N x = N y for simplification. 
Theorem 6. The difference scheme (4.6) is unconditionally stable for 1 < a, ft < 2. 
Proof. We represent the discrete functions and f™^ 1 ^ 2 into vector forms with 



77" _ („.n n n n n n . . . n n n \T 

u — ("1,1) "2,1' > U JV— 1,1) "1,2> u 2,2) ' U AT-1,2> ) "1,JV— 1: u 2,N—li ) u N-l,N-l) > 

j.n+1/2 
' JN-l,2i 



pn+1/2 _ / f n+l/2 f n+l/2 f n+l/2 ,n+l/2 f n+l/2 f n+l/2 

^ ~ Ul,l )/2,l ' ' ' ' ' /jV-1,1 ' J 1,2 ' ^2,2 ■>'''■: J ]\ 



f n+l/2 f n+l/2 f n+l/2 ,T 



and denote 



® A, + -f- 1 » At, V„ = -fs-. 



^^J^ 1 ®^*^" 1 ®^'' = W A/3 ^ / + W^ 01 ' (4 - 17) 
where the symbol <g> denotes the Kronecker product, I is the unit matrix, and matrices A a and A« are defined 
in (2.20) corresponding to a, ft, respectively. Therefore, the difference scheme (4.6) can be expressed as 

(l-V x )(l-V y )U n+1 = (I + V X )(I + V y )U n + T F n+1 l 2 , (4.18) 
then the relationship between the error e n+1 in U n+1 and the error e n in U n is given by 

e n+1 = (l-V y y\l-V x )- 1 (l + V x )(l + V y )e n . (4.19) 
Using Lemma 8, we can check that V x and V y commute, i.e., 

V x V y = V y V x = -^-^{K^Ap + K^Aj) ® (K?A a + K+A£). (4.20) 
Thus (4.19) can be rewritten as 

e n = {j y I-V y y\l + Vy)Y{{l-T) x )-\l + V x )Ye . (4.21) 
We can also calculate the symmetric part of V x by Lemma 9 as 

V x + Vl _ (K+ + K+)t i , Aq + A£ \ V y + VI _ (K£ + K^)r f A p + Aj 



2h a V 2 /' 2 2hP 

And from Theorem 2, the eigenvalues of Aa + Aa and Ap+ 2 Af) are all negative when 1 < a, ft < 2. Defining 
A a and A^ as an eigenvalue of matrices V x and V y , respectively, then it yields from the consequences of 
Lemma 4 and 7 that the real parts of A a and A^ are both less than zero. Since (1 + A a )/(1 — X a ) and 
(1 + As) /(l — Aa) are eigenvalues of matrices (I — V x )~ l (I + V x ) and (I — + respectively, 

thus the spectral radius of each matrix is less than 1, which follows that ((I — V x )~ l (I + T> x )) and 
((/ — T) y )~ x (I + Pj/)) converge to zero matrix (see Theorem 1.5 in [20]). Therefore the difference scheme 
(4.6) is unconditionally stable. □ 
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Remark 7. For the similar reason described in Remark 6 and the proof of Theorem 6, we conclude that the 
WSGD scheme with 9 weighted scheme for the time discretization for (4.1) is unconditionally stable when 



\ < e < i. 



Lemma 10. Let T> x and T> y be defined in (4.17), then 

\\{i-v x y\i-v y y l \\ 2 <i, 

IKJ-xg-^I + xgila <1, -y = x,y, 



where \\ ■ W2 denotes the 2-norm (spectral norm) 



Proof. From Theorem 2 and Lemma 7, we know that V x + T> x and V y + TfF are negative semi-definite and 



symmetric matrices. Then for any v = (v\, v%, ■ ■ ■ , v n ) T G W 1 , we obtain that 

v T v < v T (I -V^)(I -V^v, j = x,y. 
Substituting v and v T by (7 — V~ f )~ l v and v T (I — D^)" 1 , respectively, for any v G M. n , we get 

v T (I-V^y 1 (I-V y y 1 v<v T v, j = x,y. 
Thus, it leads to 

|[(7-D 7 ) 2 = sup ! — = < 1, j = x,y. 

Connsequently, 

||(/ - V x y l (I - V y )"% < ||(7 - P,)" 1 || 2 ||(7 - 2> B )-i|| 2 < 1 

holds. 

Since V x + Pj and V y + 7>J are negative semi-definite and symmetric, for any v G M. n , we have 
v T (I + P T )(7 + 7J> < v T (I - V?j){I - D>, 7 = x, y. 
By choosing vector (7 — V~ / )~ 1 v, we arrive at that for any v G M n , 

v T (7-P T )- 1 (7 + P 7 r )(7 + P 7 )(7-P 7 )- 1 « <v T t;, 7 = x,y. 

As (7 - £> 7 ) -1 (I + £7) = U + V-tW ~ ^i)" 1 ' then ^ y ields ^ 

||(7 - P 7 )~ 1 (7 + P 7 )|| 2 = ||(7 + P 7 )(7 - V^h 

; T (7 - 2?T)-i(j + 2?T)(/ + p 7 )(/ _ 2> 7 )-it> 



sup ■ „ 



< 1. 



□ 



Theorem 7. Lef it™ - Z>e exact solution of (4.1) vvz'f/i 1 < a, j3 < 2, awe? t/J^ solution of the difference 
scheme (4.6), then for all 1 < n < M, we /zctve 

||u n - C/ n || < c(r 2 + /i 2 ), (4.22) 

where c denotes the positive constant and \\ ■ \\ stands for the discrete I?-norm. 
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Proof. Let ep - = u" ■ — J7J 1 -, subtracting (4.5) from (4.6) leads to 

(/ - V x ) (I - V y )e n+1 = (I + V x ) {I + V y )e n + r£ n , (4.23) 
where V x and V y are given in 4.17 and 

e = (ei,i, e 2 ,i, • • • , 6^-1,1,61,2,62,2, • • • , ejv-1,2, ■ • • , ei.jv-i , ^2,N-\ , ' ' ' , e7v-i,Ar-i) T , 

£ = ( e l,l> £2,1, ■ ■ ■ ,£JV-1,1, £1,2,^2,2 , ' ' ' ,£JV-1,2, ' ' ' ,£l,JV-l) £2,N-1, - - - , £7V-l,AT-l) T - 

Since D a . commutes with denoting P = (I — T>^\ 1 (/ — P y ) 1 (/ + (/ + , it yields that 

e n+i = p e n + r (j _ v ^ "1 ^ _ 25 tf ) "V. (4.24) 
Iterating for all < k < n — 1 and taking the L 2 -norm on both sides, we have that 

n— 1 n— 1 

||e n || <t\\{I -V x )~\l -Vy^^^^h- ||£ n-1-fc || <r^||P fc || 2 - llf"- 1- *!, (4.25) 

fc=0 fc=0 

where Lemma 10 shows that ||(J - Da,)" 1 (I- X^)" 1 ^ < 1. 
Since and commutes, matrix P can be rewritten as 

P = (I - V^il + V X )(I - Vy)~ l (l + Vy). (4.26) 

We then obtain from Lemma 10 that 

\\P\\i < W -V X )- X {I + V X )\\ 2 \\(I -V y )- l {I + V y )\\ 2 < 1. (4.27) 
Then for any 1 < k < M, ||P fc || 2 < ||-P|| 2 < 1 holds. We can get that 

ra-l 

||e n || <r^||£ fc || <c(r 2 + /i 2 ). (4.28) 

□ 



k=0 



The convergence result for scheme (4.13) can also be obtained by the similar way as above. 

5 Numerical Examples 

5.1 One Dimensional Case 

Example 2. Consider the following problem 

= Q Dy(x,t) - e-'(:r 1+Q + r(2 + a)x), (x,t) G (0, 1) x (0, 1], (5.1) 
with the boundary conditions 

u(0,t) =0, u(l,t) = e~\ i€[0,l], 
and initial value 

u(x,0) = x l+a , x€[0,l]. 
Then the exact solution of (5.1) is u(x, t) = e~ f x 1+a . 
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Table 2: The maximum and L 2 errors and their convergence rates to Example 2 approximated by the CN- 
WSGD scheme at t = 1 for different a with r = h. 
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1.9 


16 


1.63058E-05 




2.278 14E-06 




6.02603E-05 




7.78084E-06 






32 


2.49190E-06 


2.71 


6.49029E-07 


1.81 


9.23273E-06 


2.71 


9.04790E-07 


3.10 




64 


4.93027E-07 


2.34 


1.81207E-07 


1.84 


1.35841E-06 


2.76 


1.49823E-07 


2.59 




128 


1.27340E-07 


1.95 


4.81095E-08 


1.91 


1.94436E-07 


2.80 


4.02142E-08 


1.90 




256 


3.23580E-08 


1.98 


1.24022E-08 


1.96 


3.06615E-08 


2.66 


1.12278E-08 


1.84 




512 


8.15631E-09 


1.99 


3.14892E-09 


1.98 


7.93775E-09 


1.95 


2.99226E-09 


1.91 



Example 3. Consider the following problem 

= D«u(x, t) + X D? u(x, t) + f{x, t), (x, t) G (0, 1) x (0, 1], 

u(o,t) = u(M) = o, *e[o,i], (5 - 2) 

u{x,0) = x 3 (l - xf, xe[0,l], 
with the source term 

f{x, t) = -e" 4 (* 3 (1 - x) 3 + ^ {x 3 ~ a + (1 - x) 3 " Q ) - 3 r( ^ ( ^ a) {x 4 ~ a + (1 - xf- a ) 

T(6 — a) r(7 — a) 

By simple evaluation, the exact solution of (5.2) is it(x, t) = e _t :r 3 (l — x) 3 . 

Example 4. Consider the following variable coefficients problem 

= x%D«uix, t) + (l- x) a x D?uix, t) + fix, t), ix, t) e (0, 1) x (0, 1], 

u(o,t) = u(M) = o, te[o,i], (5 - 3) 

u(x,0) = x 3 (l -x) 3 , ie[0,l], 
with the source term 

/(*,() = -e-(x 3 (l - *) 3 + ;^L(x 3 + (1 - *) 3 ) - 3j^L_(x* + (1 - x)") 

^f^y^+f'-^-^^ + a-) 6 )). 
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Table 3: The maximum and L 2 errors and their convergence rates to Example 3 approximated by the CN- 
WSGD scheme at t = 1 for different a with r = h. 









(p, q) = 


(1,0) 






(p, q) = 


(1,-1) 








II ri _ Tjn\\ 

||" u \\oo 


ftitC 


II?/™ — 77™ II 


ftitC 


\\n n — TT n II 

I" u Woo 


lcltC 


II7;™ — 7/™ll 
II" u 1 


lcltC 


1.1 


16 


1.21351E-04 


- 


6.87244E-05 


- 


1.04202E-04 


- 


5.49761E-05 


- 




32 


3.10400E-05 


1.97 


1.75798E-05 


1.97 


4.32767E-05 


1.27 


2.00595E-05 


1.45 




64 


7.93983E-06 


1.97 


4.47207E-06 


1.97 


1.48399E-05 


1.54 


7.42486E-06 


1.43 




izo 


9 01 f\lAV 06 


1 QR 

1 .70 


1 1 9QQ^P 06 


1 QR 

1 .70 


A 1 Q7RRP 06 


1 R9 
1 .oZ 


9 9^601 P 06 


1 . 1 j 




256 


5.0805 1E-07 


1.99 


2.84150E-07 


1.99 


1.10967E-06 


1.92 


6.11319E-07 


1.87 




512 


1.2751 1E-07 


1.99 


7.12580E-08 


2.00 


2.84899E-07 


1.96 


1.59692E-07 


1.94 


1.5 


16 


2.03009E-04 


- 


5.46438E-05 


- 


2.99388E-04 


- 


8.57787E-05 


- 




32 


4.52559E-05 


2.17 


1.37190E-05 


1.99 


7.90624E-05 


1.92 


2.31127E-05 


1.89 




64 


1.13225E-05 


2.00 


3.45401E-06 


1.99 


2.01483E-05 


1.97 


6.01008E-06 


1.94 




128 


2.83579E-06 


2.00 


8.67756E-07 


1.99 


5.08147E-06 


1.99 


1.53528E-06 


1.97 




256 


7.09655E-07 


2.00 


2.17555E-07 


2.00 


1.27542E-06 


1.99 


3.88274E-07 


1.98 




512 


1.77509E-07 


2.00 


5.44715E-08 


2.00 


3.19447E-07 


2.00 


9.76542E-08 


1.99 


1.9 


16 


2.02959E-04 




3.60448E-05 




2.35899E-04 




4.37067E-05 






32 


4.57927E-05 


2.15 


8.97441E-06 


2.01 


5.44882E-05 


2.11 


1.10506E-05 


1.98 




64 


9.363 12E-06 


2.29 


2.23928E-06 


2.00 


1.13848E-05 


2.26 


2.77301E-06 


1.99 




128 


2.03859E-06 


2.20 


5.59714E-07 


2.00 


2.55286E-06 


2.16 


6.94607E-07 


2.00 




256 


5.08948E-07 


2.00 


1.39944E-07 


2.00 


6.35234E-07 


2.01 


1.73827E-07 


2.00 




512 


1.27160E-07 


2.00 


3.49898E-08 


2.00 


1.58420E-07 


2.00 


4.34792E-08 


2.00 



By simple evaluation, the exact solution of (5.3) is u(x, t) = e x (1 — x) . 



5.2 Two Dimensional Case 

Example 5. The following fractional diffusion problem 

du(x,y,t) 



Of 



Q D l £ 2 u(x, y, t) + x D\ 2 u{x, y, t) + Dl s u(x, y, t) + y D\ % u(x, y, t) + f(x, y, t) 



is considered in the domain Q = (0, l) 2 and t > with boundary conditions u(x,y,t)\on = and the 
initial condition u(x, y, 0) = x 3 (l — x) 3 y 3 (l — y) 3 , where the source term 



f(x,y,t) 



(x 3 (i - x)V(i - yf) + (M-i* 1 - 8 + (i - *n - §^(x 2 - 8 + (i - -) 2 - 8 ) 



+ 



+ 



r(4.i 
r(4) 
r(2.2) 



+ (l-x) 3 - 8 ) 



r(2. 
r(7) 



r(3.8) 

+ (l-x) A - s ))y 3 (l-y) 3 



V 2 +u-y) 3 - 2 ) 



r(5. 

-^(^ + d-y)") 

r(7) , 4 . 2 



^ + (l-y) 4 ' 2 ))x 3 (l-x) 



r(4.2) w v ^ y r(5.2) 

r/ie« etacf solution of the fractional partial differential equation is u(x,y,t) = e~ t x 3 (l — x) 3 y 3 (l — y) 3 . 

We use four numerical schemes: LOD (4.14), PR-ADI (4.7), Douglas- ADI (4.8) and D'yakonov-ADI 
(4.9), to simulate Example 5, the maximum and L 2 errors and their convergence rates to Example 5 ap- 
proximated at t = 1 are listed in Table 5, where N = N x = N y , and p, g are the shifted numbers of the 
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Table 4: The maximum and L 2 errors and their convergence rates to Example 4 approximated by the CN- 
WSGD scheme at t = 1 for different a with r = h. 









(p, q) = 


(1,0) 






(p, q) = 


(1,-1) 






JV 


II n _ Tjn\\ 
II" u 1 1 oo 


I"£ltC 


II?/" — T7 n \\ 
II a u | 




\\n n — TT n \\ 

I" u Woo 


1'B.tC 


II" u 1 


TcltC 


1.1 


16 


1.77123E-04 




7.32001E-05 




3.95613E-04 




1.92219E-04 










1 .70 


1 7A1 RzLP 


1 fK 


y. / J /OjlL-KJJ 


L.\jL 




L.LL 




64 


1.08962E-05 


2.04 


4.36356E-06 


2.01 


2.43654E-05 


2.00 


1.00363E-05 


2.04 




128 


2.66784E-06 


2.03 


1.08906E-06 


2.00 


6.10991E-06 


2.00 


2.51523E-06 


2.00 




256 


6.67126E-07 


2.00 


2.72235E-07 


2.00 


1.53026E-06 


2.00 


6.31764E-07 


1.99 


1.5 


16 


1.88510E-04 


- 


6.18902E-05 


- 


3.56874E-04 


- 


1.30433E-04 


- 




32 


4.48741E-05 


2.07 


1.46628E-05 


2.08 


8.32954E-05 


2.10 


2.80619E-05 


2.22 




64 


1.10524E-05 


2.02 


3.61334E-06 


2.02 


2.02076E-05 


2.04 


6.65178E-06 


2.08 




128 


2.74933E-06 


2.01 


8.99424E-07 


2.01 


4.98975E-06 


2.02 


1.63398E-06 


2.03 




256 


6.86120E-07 


2.00 


2.245 18E-07 


2.00 


1.24092E-06 


2.01 


4.05976E-07 


2.01 


1.9 


16 


1.61881E-04 




4.02897E-05 




1.79407E-04 




5.75044E-05 






32 


3.43080E-05 


2.24 


9.582 13E-06 


2.07 


4.06728E-05 


2.14 


1.27751E-05 


2.17 




64 


7.72475E-06 


2.15 


2.35289E-06 


2.03 


9.30708E-06 


2.13 


3.02268E-06 


2.08 




128 


1.91676E-06 


2.01 


5.83977E-07 


2.01 


2.34573E-06 


1.99 


7.37420E-07 


2.04 




256 


4.80573E-07 


2.00 


1.45527E-07 


2.00 


5.9261 1E-07 


1.98 


1.82315E-07 


2.02 



WSGD operators. From the numerical results, three ADI schemes obtain more accurate solution than the 
LOD scheme, and it also reflects that the three ADI schemes are equivalent in two dimensional case. 

6 Conclusion 

The paper provides the novel second order approximations for fractional derivatives, called the weighted and 
shifted Griinwald difference operator; it also suggests a direction to gain higher order discretization and com- 
pact schemes of fractional derivatives. The discretizations are used to solve one and two dimensional space 
fractional diffusion equations; several numerical schemes are designed, their effectiveness are theoretically 
proved and numerically verified. 
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